Resolving the HBT Puzzle in Relativistic Heavy Ion Collisions 
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Two particle correlation data from the Relativistic Heavy Ion Collider have provided detailed 
femtoscopic information describing the space-time structure of the emission of pions. This data had 
avoided description with hydrodynamic-based approaches, in contrast to the success of hydrodynam- 
ics in reproducing other classes of observables. This failure has inspired the term "HBT puzzle", 
where HBT refers to femtoscopic studies which were originally based on Hanbury-Brown Twiss in- 
terferometry. Here, the puzzle is shown to originate not from a single shortcoming of hydrodynamic 
models, but the combination of several effects: mainly including pre-thermalized acceleration, using 
a stiffer equation of state, and adding viscous corrections. 



Experiments at the Relativistic Heavy Ion Collider 
(RHIC) have revealed a new state of matter, the 
strongly interacting quark gluon plasma, which ap- 
pears to have both perhaps the lowest ratio of vis- 
cosity to entropy of any measured substance [l], Q. 
This conclusions is based on comparisons of hydro- 
dynamic models with experimental spectra and large 
angle correlations that reveal strong radial and ellip- 
tic collective flow, i.e., flow relative to the original 
beam axis. Whereas spectra and large angle correla- 
tions are consistent with ideal hydrodynamics [H, 0] , 
hydrodynamic models have poorly reproduced corre- 
lations at small-relative momentum [4|, [5|, l6|. These 
correlations are related to the spatial and temporal 
properties of pion emission 0], and are often referred 
to a HBT (Hanbury-Brown Twiss) measurements af- 
ter similar measurements with light 0. It appears 
that hydrodynamic models underestimate the explo- 
siveness of the collision, or equivalently overestimate 
the duration of the emission process. In contrast, 
some microscopic approaches have been more success- 
ful, though not completely satisfactory, in reproduc- 
ing the data 0, 0, El [H| • Unlike the hydrodynamic 
models, which employed first order phase transitions, 
the effective equations of state for the microscopic ap- 
proaches are extremely stiff, leading to more explosive 
collisions. In short, the HBT puzzle involves finding 
whether one can reproduce femtoscopic observations 
with hydrodynamic models without using an equation 
of state that is inconsistent with lattice QCD calcu- 
lations. In this study, we show this can be accom- 
plished if three improvements are incorporated into 
hydrodynamic models: accounting for the buildup of 
collective flow in the first instants of the collision be- 
fore thermalization is attained, using a stiffer equation 
of state, and including viscosity. The hydrodynamic 
model used here is outlined in detail in Ref . [l3| . 



The Koonin equation [14j relates the experimentally 
measured correlation function to the outgoing phase 



space density, 

C(P,q)= / d\ S(P,r)|<Xq,r)| 2 . (1) 

Here, P and q are the total and relative momentum. 
The source function S(P,r) describes the probability 
for two particles with identical momenta k = P/2, 
to be separated by r in their asymptotic state if the 
relative interaction between the particles were to be 
ignored. Any dynamical model, whether based on hy- 
drodynamics or microscopic degrees of freedom, pro- 
vides a list of positions and times from which particles 
of specific k are emitted, and can then be used to gen- 
erate the source function. 

With some effort, the source function can be imaged 
from the correlation function, but for the purposes of 
this paper we will only consider parameters extracted 
by fitting to correlations that arise from a Gaussian 
source, 



S(P,r) ~exp 
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The parameter i?i on g describes the longitudinal size 
along the beam axis, i? s ide describes the extent along 
the sideward dimension perpendicular to both the 
beam axis and P. The outward dimension, i? ut de- 
scribes the size along the axis along the direction of 
P, assuming one has boosted along the beam axis to 
a frame where P points perpendicular to the x axis. 
Each of these dimensions can be a function of the 
transverse momentum k t = P x /2 or the longitudinal 
rapidity of the pair. For central collisions, which is 
what will be described here, the dimensions are inde- 
pendent of the azimuthal angle of P. Furthermore, we 
will only consider central rapidity and only consider 
the kt dependence of the Gaussian dimensions. 

Figure [1] shows experimentally determined radius 
parameters from 100 A GeV Au on 100 A GeV Au col- 
lisions at RHIC. For comparison, source sizes were 
generated from a hydrodynamic model coupled to a 
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FIG. 1: (color online) Gaussian radii reflecting spatial sizes 
of outgoing phase space distributions in three directions: 
Rout, Rside and -Riong- Data from the STAR collaboration 
(red stars) are poorly fit by a model with a first-order 
phase transition, no pre-thermal flow, and no viscosity 
(solid black squares). Correcting for all those deficiencies, 
and using a more appropriate treatment of the relative 
wave function in Eq. (JXJ brings calculations close to the 
data (filled black circles). The sequential effects of includ- 
ing prethermal acceleration (open blue squares), using a 
more realistic equation of state (open green diamonds), 
and adding viscosity (open cyan triangles) all make sub- 
stantial improvements to fitting the data. An improved 
relative wave function yielded modest improvements (com- 
pare open cyan triangles to filled black circles). 



cascade code. The cascade microscopically simulates 
the final stages of the collision and breakup where lo- 
cal kinetic equilibrium is lost and hydrodynamics is 
unjustified. The times and positions of last collisions 
for particles of a specific k were used to calculate the 
source function, from which correlation functions were 
generated via Eq. |T]). These were then fit to corre- 
lations from Gaussian sources to extract radii, which 
are also displayed in Fig. [1] 

As a benchmark, the first calculation (filled squares 
in Fig. [T]) was parameterized similarly to previous hy- 
drodynamic calculations, and failed in a similar man- 
ner. Transverse expansion was delayed until 1 fm/c 



after the initial collision. A strong first-order phase 
transition, which is inconsistent with lattice gauge 
theory, was employed, and the viscosities were set to 
zero. Additionally, an over-simplified relative wave 
function, neglecting Coulomb and strong interactions 
between the pions, was used to generate correlation 
functions. Since the source functions are not truly 
Gaussian, this can lead to different Gaussian radii. 
This benchmark calculation overstates the R out / i? s idc 
ratio by ~ 40% and overstates i?i ong by ~ 25%. 

The second calculation (open squares in Fig. [T|) ac- 
counts for prethermal acceleration by beginning the 
expansion 0.1 fm/c after the initial collision, roughly 
the amount of time required for the Lorentz contracted 
nuclei to traverse one another. The importance of 
pre-thermalized acceleration has been emphasized in 
several studies during the last few years [l(J 15. 
As was shown in Ref. 117 1. flow during the first 1 



fm/c is approximately universal for any system with 
a traceless energy tensor, including partonic and field 
based pictures, independent of thermalization. Since 
the transverse expansion starts earlier, the longitudi- 
nal size is smaller at breakup, more in line with data. 
The Rout /Rside ratios also drop, moving modestly to- 
ward the data. 

The second improvement to be considered is to use 
a stiffer equation of state. Early studies used an equa- 
tion of state with a first order phase transition with 
a large latent heat [1, [H, 0|. Such soft equations of 
state have constant temperature and pressure for en- 
ergy densities between eh and e/,. + L, where eh is the 
maximum density of the hadronic phase. Here, eh 
corresponds to a hadronic gas with a temperature of 
T c = 170 MeV, and L is the latent heat. In con- 
trast, lattice QCD now suggests a crossover transition 
where the pressure rises continuously with energy den- 
sity. There indeed exists a soft region, but the speed 
of sound, c 2 s = dP/de, never falls below 0.1 and the 
width of the soft region is somewhat lower than the 
latent heat L assumed in the previous studies. The 
benchmark calculation, displayed in the upper panel, 
assumed a first order transition with a latent heat 
L = 1.6 GeV/fm 3 with a lower bound to the mixed 
phase at eh ~ 500 MeV/fm 3 . This is not only inconsis- 
tent with lattice calculations, but is also inconsistent 
with femtoscopic analyses of data at lower energies. 
For heavy ion collisions at the upper AGS and for the 
lower SPS beam energies, maximum energy densities 
were in the neighborhood of eh + L. For a first order 
phase transition the pressure P stays fixed through- 
out the mixed phase, and these conditions would have 
minimal values of P/e with minimal explosivity re- 
sulting in perhaps dramatically large lifetimes, well 
exceeding 20 fm/c. The long duration of the emis- 
sion would lead to extended values of the outward 
dimensions of the phase-space cloud [HI, 0|. This 
was not observed. The third calculation (open dia- 
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monds in Fig. [T]) assumed a soft region of half the 
width in energy density, and with a speed of sound 
of Cj = 0.1, rather than zero for a first-order transi- 
tion. Once above the soft region, both calculations 
assumed a stiffening with the speed of sound c? s = 0.3. 
One can consult Ref. [2(| for a more sophisticated at- 
tempt at parameterizing lattice equations of state. As 
expected, the stiffer equation of state leads to more 
explosive collisions. A more sudden emission results 
in phase space clouds, f(k t ,r,t — » oo), that are less 
extended along the outward dimension, which lowers 
the the -R ut/-Rsidc ratio. As can be seen in Fig. [1] 
the ratio moves significantly toward the data with the 
stiffer equation of state. 

Shear viscosity is also known to increase the ex- 
plosiveness of the collision 21, This can be 



understood by considering viscous corrections to the 
stress-energy tensor. At early times the velocity gra- 
dient is largely longitudinal, which affects the stress- 
energy tensor by strengthening the transverse pres- 
sure, T xx = T yy , and decreasing the longitudinal pres- 
sure. In the Navier-Stokes equation, 
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where rj is the shear viscosity and the velocity gradi- 
ent for early times is dv z /dz = 1/t. After the first 
few fm/c, the transverse acceleration is determined 
by T xx and T yy . As originally demonstrated in |22j, 
and shown here in the lower panel of Fig. [TJ the 
^out/-Rsidc ratio can be lowered by ~ 10% with re- 
alistic shear viscosities. Analyses of elliptic flow have 
pointed to a small shear viscosity, perhaps approach- 
ing the KSS limit [231 ]. 7y m i n = s/Att, where s is the 
entropy density. The neglect of pre-thermalized flow 
in these calculations might have led to underestimates 
of the viscosity, but nonetheless, it is expected that 77 
is not much greater than the KSS bound. Below T c , 
collisions are binary and the cascade prescription nat- 
urally accounts for viscous effects. Bulk viscosity is 
expected to be important near T c due to the inabil- 
ity of the system to maintain equilibrated chiral fields 
near T c 24]. These expectations were also verified 
with lattice calculations [25j]. The impact of adding 
viscosity is shown in Fig. [1] (open triangles). Shear, 
which was set at twice the KSS bound, significantly 
affected the sources sizes, while changes to the bulk 
viscosity had little effect. Combined with the previ- 
ously discussed effects, the calculation now approach 
the experimental data. 

As a final improvement to the analysis, a more real- 
istic scattering wave function was used in Eq. (p} . The 
resulting correlation function was then fit to a Gaus- 
sian source using the approximate treatment applied 
in experimental treatments [2^, |27j|. Since sources 
are not purely Gaussian, and since the experimen- 
tal fitting procedure only approximates the effects of 



Coulomb in the wave functions, the resulting radii dif- 
fer slightly. The final results (filled circles in Fig. [1]) 
now match the data within rj 5%, which is less than 
the normally quoted systematic experimental errors 
of 5-10%. Thus, the HBT puzzle appears to have re- 
sulted from a conspiracy between several shortcomings 
in the original models, each of which led to an under- 
estimate of the cxplosiveness of the collision, and each 
of which pushed i? ut/-Rsidc in the same direction. 

Other adjustments to the dynamics are known to af- 
fect femtoscopic radii. For the calculations presented 
here, the initial energy density profile was scaled to 
fit experimental charge multiplicity, while the shape of 
the profile was set by the wounded nucleon model [28| . 
More compact profiles have been shown to increase the 
explosiveness and lower the i?out/-R S idc ratio 13. 
The model used here assumed boost-invariant acceler- 
ationless flow along the beam axis, similar to a simple 
Hubble expansion. This symmetry simplified and ac- 
celerated hydrodynamic calculations, but is only ap- 
proximate, and is expected to cause errors of a few 
percent in HBT radii, most likely leading to a lower- 
ing of i?ion g by a few percent [29(. Another effect not 
considered here concerns the interaction of the outgo- 
ing pions with a mean field from the remaining matter 
[30(. Such effects only come into play after the final 
collisions of the pions, when densities are < 0.1 per 
fm 3 and fields are expected to be small. Mean field 
distortions are expected to affect radii only at the level 
of a few percent 31 1 . 

One of the puzzling aspects of HBT analyses came 
from fits with bast-wave models, which are parame- 
teric pictures based on thermal emission from a col- 
lectively expanding source. Parameters include a 
breakup temperature T, an outer radius R, a breakup 
time r, a linearly rising transverse collective veloc- 
ity with a maximum, v ma x, and an emission duration 
At. Blast wave fits are similar in quality to what was 
achieved here with a dynamic model. The surprising 
outcome of the fits was that radii were on the order of a 
dozen fm, with emission confined to within a few fm/c 
of t = 10 fm/c [32j, 133J. These parameters suggest an 
unphysically high breakup density. Figure[2]shows the 
outward coordinate x and the time r at which emis- 
sion occured for the final model in Fig. [1] and are 
similar to what was seen in 12j. Points are shown 



only for particles emitted with momentum p x = 300 
MeV/c. The dashed lines have a slope correspond- 
ing to the velocity v x = p x / E. Any two emissions at 
points connected by such a line contribute identically 
to the final phase space density, and thus identically 
to the femtoscopy. The positive correlation between x 
and t allows longer-lived emission to femtoscopically 
mimic more sudden emission at a much earlier time. 
This emphasizes the importance of using realistic dy- 
namical models to compare to femtoscopic data and 
underscores the limits of parametric fits. 



4 




-20 



20 

x (fm) 



40 



FIG. 2: (color online) Final emission positions and times 
for particles with transverse momentum of 300 MeV/c 
along the x axis. Emission has a positive correlation be- 
tween position and time, though lags behind the slope of 
the velocity (illustrated by dashed lines). Due to the posi- 
tive x — t correlation, emissions of longer duration can still 
result in phase space clouds that are compact along the 
outward direction, with -R ou t/-Rsidc ~ 1. 



The calculations shown here were successful in 
matching experimental values of the mean transverse 
momentum for pions, kaons and protons. Elliptic flow, 
analyses of which are inherently precluded with the 
model used here due to an assumption of azimuthal 
symmetry, needs to be fit with the same parameters. 
More detailed aspects of femtoscopic data also need 
to be matched such as: radii with respect to the reac- 
tion plane [2I], correlations of other species, especially 
non-identical particles 3^, 36 L and non-Gaussian fea- 
tures of the source function [37J . Even if all these data 
are reproduced, it does not fully validate the model. 
That would require an ambitious statistical analysis 
of the set of model parameters and assumptions, sim- 
381 ] . Although these goals require significant 



ilar to 

effort in the coming years, the current analysis has 
eliminated any puzzle about femtoscopy for the time 
being, as the experimental radii appear to be satisfac- 
torily described within a rather standard theoretical 
picture of how RHIC collisions evolve. 
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